rm(list=ls())
library(fields); library(maps); library(ncdf4);library(abind)
library(geosphere);library(mapdata)
setwd('/Volumes/Disk2/JJA_ozone_MAM/Harvard_Dataverse')

#----- load functions -----------------------
source("./Function/get_geo.R")
source("./Function/get_met.R")
source('./Function/read_detrend.R')
source('./Function/eof.mca.R')
source('./Function/cov4gappy.R')
source('./Function/Henderson.R')
source('./Function/filled.contour3.R')

#-----------------------------------------
ss1=load("./Figures/Figure 2/mov_detrend.Rdata")

col.gen = colorRampPalette(c("#0563FA", "white","#ff4d4d"), space="rgb")
levels=seq(-0.7, 0.7, length=31)
mid=(levels[1:(length(levels)-1)]+levels[2:(length(levels))])/2
cols=col.gen(length(levels)-1)
if(min(abs(mid))>0){
	ind1=which(mid<=0);ind2=which(mid>=0)
	ind=abind(Reduce(intersect, list(ind1+1, ind2)), Reduce(intersect, list(ind1, ind2-1)))
	ind=sort(unique(ind))
  	 cols[ind]="white"
}	

#======= Part 1: read data =======================================
mai=c(0.1, 0.05, 0.2, 0.05)
mgp=c(1.2, 0.4, 0)
tcl=-0.2
ps=12

close.screen(all.screens = TRUE)
m <- rbind(c(0, 0.33, 0.5, 0.9), c(0.33, 0.67, 0.5, 0.9), c(0.67, 1.0, 0.5, 0.9), c(0.2,0.8,0.43,0.68))
split.screen(m)

screen(1)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
result=find.cor3(yearly_O3, PC1)
ap=result$corr
p=result$p
ap[ap>=0.7]=0.7;ap[ap<=-0.7]=-0.7
image(lon.frame[14:26],lat.frame, ap[14:26,],zlim=c(-0.7,0.7),col=cols,xaxt='n',yaxt='n',xlab='',ylab='');
map("world", add = TRUE)
cr=0.05;p[is.na(p)]=1
 for (i in 1:28){
	 for (j in 1:11){
		 if(abs(p[i,j])<cr)points(lon.frame[i],lat.frame[j],pch=16,cex=0.4,col=1)
	 }
 }
 title(expression(paste('(a) Corr ozone vs. MAM-',Delta,'SST' ,sep='')),cex.main=1.1)
box()

screen(2)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
result=find.cor3(yearly_O3, PC2)
ap=-result$corr
p=result$p
ap[ap>=0.7]=0.7;ap[ap<=-0.7]=-0.7
image(lon.frame[14:26],lat.frame, ap[14:26,],zlim=c(-0.7,0.7),col= cols,xaxt='n',yaxt='n',xlab='',ylab='');
map("world", add = T)
 cr=0.05;p[is.na(p)]=1
 for (i in 1:28){
	 for (j in 1:11){
		 if(abs(p[i,j])< cr)points(lon.frame[i],lat.frame[j],pch=16,cex=0.4,col=1)
	 }
 }
 title(expression(paste('(b) Corr ozone vs. MAM-',Delta,'SLP' ,sep='')),cex.main=1.1)
box()

screen(3)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
	ap=array(NA,c(dim(yearly_O3)[1],dim(yearly_O3)[2]))	
	p_values=array(NA,c(dim(yearly_O3)[1],dim(yearly_O3)[2]))	
	for(ii in 1:dim(yearly_O3)[1]){
		for(jj in 1:dim(yearly_O3)[2]){
			if(sum(!is.na(predict_O3[ii,jj,]))>=5){
			fit<-cor.test(yearly_O3[ii,jj,], predict_O3[ii,jj,])
			ap[ii,jj]=fit$estimate
			p_values[ii,jj]=fit$p.value
			}
		}
	}
ap[ap>=0.7]=0.7;ap[ap<=-0.7]=-0.7
image(lon.frame[14:26],lat.frame, ap[14:26,],zlim=c(-0.7,0.7),col= cols,xaxt='n',yaxt='n',xlab='',ylab='');
map("world", add = T)
cr=0.05; p_values[is.na(p_values)]=1
 for (i in 1:28){
	 for (j in 1:11){
		 if((p_values[i,j])<cr & ap[i,j]>0)points(lon.frame[i],lat.frame[j],pch=16,cex=0.4,col=1)
	 }
 }
title('(c) Cross-validated Correlation',cex.main=1.1,font.main=1)
box()

screen(4)
mai=c(0.1, 0.1, 0.05, 0.2)
par(mai=mai, mgp=mgp, tcl=tcl, ps=12)
image.plot(ap,legend.shrink=0.9,legend.only=TRUE,zlim=c(-0.7,0.7),col= cols,horizontal=TRUE,legend.width=2.8,axis.args = list(mgp = c(0, 0.8, 0),tcl=-0.2,padj=-1))
text("r",x=par("usr")[2]-0.23,y=par("usr")[3]+0.485,srt=0, adj = 0,xpd=TRUE,cex=1.5)
